close all; clear all; clc;
mylabels = {'H. rufescens', ... % 1
            'L. variegatus', 'S. purpuratus', 'A. punctulata', 'L. pictus', ... % 2:5
            'C. intestinalis', ... % 6
            'D. rerio', ... % 7
            'B. taurus', ... % 8 
            'H. sapiens', 'H. sapiens (viscous)', ... % 9:10
            };
mycolors = [0 1 1;
            1 0 0;
            1 0 0;
            1 0 0;
            1 0 0;
            1 0 1;
            0 1 0;
            0 0 1;
            0 0 1;
            0 0 1];  
mysymbol = 'os^vooos^vdpo';
mymarkersize = 3;
myfontsize = 9;


%% load data
% distance matrix
load('200617KDmx-SpecToSpec.mat')
    distancematrix = wMij;
% flip l. var and s. purp
    distancematrix(2:3,:) = [distancematrix(3,:);distancematrix(2,:)];
    distancematrix(:,2:3) = [distancematrix(:,3) distancematrix(:,2)];
% flip a punc and l pic
    distancematrix(4:5,:) = [distancematrix(5,:);distancematrix(4,:)];
    distancematrix(:,4:5) = [distancematrix(:,5) distancematrix(:,4)];
        
%% FIGURE 3
%% Figure 3a: plot distance matrix
ind_keep = 1:9;
    distancematrix_now  = distancematrix(ind_keep,ind_keep);
    mylabels_now = mylabels(ind_keep);
    distancematrix_minmax = [0 max(distancematrix_now(:))];
ind_order = [1 5 2 3 4 6 7 8 9];
    distancematrix_now = distancematrix_now(ind_order,:);
    distancematrix_now = distancematrix_now(:,ind_order);
    mylabels_now = mylabels_now(ind_order);
mypapersize = [8.9 8];
    AR = 1; h_off = 0.2*5; w_off = 2.4;  h = 5; w = h*AR; d = 0.1;
h_fig = figure;
    set(h_fig,'Units', 'centimeters', 'Position', [1 1 mypapersize],'PaperUnits', 'centimeters', 'papersize', mypapersize,'Color','w');
    ax_distancematrix = axes('Units','centimeters','Position', [w_off, h_off, w, h]);
imagesc(distancematrix_now, distancematrix_minmax); colormap('hot');
    set(ax_distancematrix, 'xaxislocation', 'top');
    set(ax_distancematrix, 'XTick', 1:length(ind_keep), 'YTick', 1:length(ind_keep));
    set(ax_distancematrix, 'XTickLabel', mylabels_now, 'fontangle', 'italic', 'fontsize', myfontsize, ...
                           'YTickLabel', mylabels_now, 'fontangle', 'italic', 'fontsize', myfontsize);
    xtickangle(ax_distancematrix, 90);
    ax_cbar = axes('Units','centimeters','Position', [w_off, h_off-d-0.25, w, 0.25]);
    imagesc(linspace(distancematrix_minmax(1),distancematrix_minmax(2),128),'XData',linspace(distancematrix_minmax(1),distancematrix_minmax(2),128),'YData',zeros(1,100));axis xy;
    set(gca, 'YTick', [], 'fontsize', myfontsize);
    set(ax_cbar);%, 'xaxislocation', 'bottom')
    xlabel('Kinematic distance, d', 'fontsize', myfontsize);
  

%% Figure 3b: plot kinematic tree
mypapersize = [10 7];
    d = 0.2; AR = 4/3;h_off = 1.;w_off = 1.; h = mypapersize(2)-h_off-d; w = h*AR; 
ind_keep = 1:9;
    mynames= mylabels(ind_keep);
    distancematrix = distancematrix([ind_keep],[ind_keep]);
kinematictree = seqlinkage(distancematrix,'single',mylabels(ind_keep));
    kinematictree = reorder(kinematictree, [9 7 8 5 6 3 4 1 2]);
plot(kinematictree,'Orientation','left')
h_fig = gcf;
  set(h_fig,'Color','w');
  set(h_fig,'Units', 'centimeters', 'Position', [1 1 mypapersize],'PaperUnits', 'centimeters', 'papersize', mypapersize);
h_tree = get(h_fig,'userdata');
    set(h_tree.BranchDots,'Marker', 'none');
    set(h_tree.LeafDots,'Marker', 'none'); 
    set(h_tree.BranchLines,'LineWidth', 1.3); 
    set(h_tree.axes,'Units','centimeters');
        pos = get(h_tree.axes,'position');
        pos = [0.4611    0.6504    6    5.8366];
        set(h_tree.axes,'position', pos);    
axis off


%% Figure 3c: plot full 18s tree
load('matrixUPGMA.mat');
my18stree = seqlinkage(distanceMatrix_18s,'single',names_18s);
    my18stree = reorder(my18stree, [9 8 6 5 7 4 3 2 1]);
plot(my18stree,'Orientation','right');
h_fig = gcf;
  set(h_fig,'Color','w')
  set(h_fig,'Units', 'centimeters', 'Position', [1 1 mypapersize],'PaperUnits', 'centimeters', 'papersize', mypapersize);
h_tree = get(h_fig,'userdata');
    set(h_tree.BranchDots,'Marker', 'none');
    set(h_tree.LeafDots,'Marker', 'none');
    set(h_tree.BranchLines,'LineWidth', 1.3);
    set(h_tree.axes,'Units','centimeters');
        pos = get(h_tree.axes,'position');
        pos = [pos(1)+pos(3)-6    0.6504    6    5.8366];
        set(h_tree.axes,'position', pos)    ;
axis off;
   